/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __NOMINAL_HH__
#define __NOMINAL_HH__

#include <cstring>
#include <vector>

#include "fastonebigheader.hh"
#include "model.hh"
#include "rating.hh"

namespace playerpiano
{
  class Nominal
    {
      private:
        unsigned int n_workers;
        unsigned int n_labels;
        float* alpha;
        float* gamma;
        float prioralpha_sq_stddev;

        static inline float
        square (float x)
          {
            return x * x;
          }

        struct Alphafunc
          {
            const float* alpha;
            unsigned int true_label;

            Alphafunc (const float* _alpha,
                       unsigned int _true_label) 
              : alpha (_alpha),
                true_label (_true_label)
              {
              }

            inline float
            operator() (unsigned int k) const
              {
                return (k == true_label) ? -1.0 : alpha[k];
              }
          };

        void
        dlogpratingdalpha (float*        dalpha,
                           const Rating& rating,
                           unsigned int  true_label,
                           float         beta)
          {
            Alphafunc alphafunc (alpha + rating.worker * n_labels * n_labels
                                       + true_label * n_labels,
                                 true_label);
            unsigned int maxindex = 0;

            for (unsigned int i = 1; i < n_labels; ++i)
              {
                if (-alphafunc (i) > -alphafunc (maxindex) )
                  {
                    maxindex = i;
                  }
              }

            float max = -alphafunc (maxindex) * beta;
            float denominator = 0;
            float ddenominator_dalpha[n_labels];

            memset (ddenominator_dalpha,
                    0,
                    n_labels * sizeof (ddenominator_dalpha[0]));

            for (unsigned int i = 0; i < n_labels; ++i)
              {
                float d = -alphafunc (i) * beta - max;
                float expd = fastexp (d);
                denominator += expd;
                ddenominator_dalpha[i] -= expd * beta;
                ddenominator_dalpha[maxindex] += expd * beta;
              }

            for (unsigned int i = 0; i < n_labels; ++i)
              {
                dalpha[i] = -ddenominator_dalpha[i] / denominator;
              }

            dalpha[rating.label] -= beta;
            dalpha[maxindex] += beta;
          }

      public:
        Nominal (float*       _alpha,
                 float*       _gamma,
                 float        prioralpha_stddev,
                 unsigned int _n_workers,
                 unsigned int _n_labels)
          : n_workers (_n_workers),
            n_labels (_n_labels),
            alpha (_alpha),
            gamma (_gamma),
            prioralpha_sq_stddev (square (prioralpha_stddev))
          {
          }

        Nominal (Model* model)
          : n_workers (1U << strtoul (static_cast<char*> (model->lookup ("--n_worker_bits").first), 0, 0)),
            n_labels (strtoul (static_cast<char*> (model->lookup ("--n_labels").first), 0, 0)),
            alpha (static_cast<float*> (model->lookup ("alpha").first)),
            gamma (static_cast<float*> (model->lookup ("gamma").first)),
            prioralpha_sq_stddev (1.0)
          {
          }

        void
        apply_deferred_prior (unsigned int  w,
                              uint64_t      s,
                              uint64_t      t,
                              float         eta,
                              float         rho,
                              unsigned int  n_items)
          {
            if (t > s)
              {
                float powt = fastpow (t, 1.0f - rho);
                float sigma = prioralpha_sq_stddev; 
                float pows = fastpow (s, 1.0f - rho);
                float decay = fastexp (eta * (pows - powt) / 
                                       (n_items * (1.0f - rho) * sigma));

                float* alphaptr = alpha + w * n_labels * n_labels;

                for (unsigned int j = 0; j < n_labels; ++j)
                  {
                    for (unsigned int k = 0; k < n_labels; ++k)
                      {
                        if (k == j)
                          {
                            /* do nothing */
                          }
                        else
                          {
                            float dalphatmp = 
                              (1.0 - decay) * 
                              (gamma[j * n_labels + k] - 
                               alphaptr[j * n_labels + k]);

                            alphaptr[j * n_labels + k] += dalphatmp;
                            gamma[j * n_labels + k] += 
                              dalphatmp / (1.0 + n_workers);
                          }
                      }
                  }
              }
          }

        void
        update (uint64_t                  t,
                float                     eta,
                float                     rho,
                float                     beta,
                const std::vector<float>& pz,
                const Rating&             rating)
          {
            float* alphaptr = alpha + rating.worker * n_labels * n_labels;
            float etat = eta * fastpow (t, -rho);

            for (unsigned int j = 0; j < n_labels; ++j)
              {
                float tmpdalpha[n_labels];

                dlogpratingdalpha (tmpdalpha, rating, j, beta);

                for (unsigned int k = 0; k < n_labels; ++k)
                  {
                    if (k == j)
                      {
                        /* do nothing */
                      }
                    else
                      {
                        alphaptr[j * n_labels + k] += 
                          etat * pz[j] * tmpdalpha[k];
                        gamma[j * n_labels + k] +=
                          etat * pz[j] * tmpdalpha[k] / (1.0 + n_workers);
                      }
                  }
              }
          }

        float
        logprating (const Rating&   rating,
                    unsigned int    true_label,
                    float           beta) const
          {
            Alphafunc alphafunc (alpha + rating.worker * n_labels * n_labels
                                       + true_label * n_labels,
                                 true_label);
            unsigned int maxindex = 0;

            for (unsigned int i = 1; i < n_labels; ++i)
              {
                if ( -alphafunc (i) > -alphafunc (maxindex) )
                  {
                    maxindex = i;
                  }
              }

            float max = -alphafunc (maxindex) * beta;
            float denominator = 0;

            for (unsigned int i = 0; i < n_labels; ++i)
              {
                float d = -alphafunc (i) * beta - max;
                float expd = fastexp (d);
                denominator += expd;
              }

            return -alphafunc (rating.label) * beta - max - fastlog (denominator);
          }

        float
        logprating (const Rating&   rating,
                    unsigned int    true_label,
                    float           beta,
                    float&          dlogbeta,
                    float&          ddlogbeta) const
          {
            Alphafunc alphafunc (alpha + rating.worker * n_labels * n_labels
                                       + true_label * n_labels,
                                 true_label);

            unsigned int maxindex = 0;

            for (unsigned int i = 1; i < n_labels; ++i)
              {
                if ( -alphafunc (i) > -alphafunc (maxindex) )
                  {
                    maxindex = i;
                  }
              }

            float max = -alphafunc (maxindex) * beta;
            float denominator = 0;
            float ddenominatordlogbeta = 0;
            float dddenominatorddlogbeta = 0;

            for (unsigned int i = 0; i < n_labels; ++i)
              {
                float d = -alphafunc (i) * beta - max;
                float expd = fastexp (d);
                denominator += expd;
                ddenominatordlogbeta += expd * d;
                dddenominatorddlogbeta += expd * d * (d + 1.0);
              }

            dlogbeta = -alphafunc (rating.label) * beta - max
                        - ddenominatordlogbeta / denominator;
            ddlogbeta = -alphafunc (rating.label) * beta - max
                        + square (ddenominatordlogbeta / denominator)
                        - dddenominatorddlogbeta / denominator;

            return -alphafunc (rating.label) * beta - max - fastlog (denominator);
          }

        static void
        ensure (Model* model)
          {
            unsigned int n_workers (1U << strtoul (static_cast<char*> (model->lookup ("--n_worker_bits").first), 0, 0));
            unsigned int n_labels (strtoul (static_cast<char*> (model->lookup ("--n_labels").first), 0, 0));
            model->ensure ("alpha",
                           n_workers * square (n_labels) * sizeof (float));
            model->ensure ("gamma", square (n_labels) * sizeof (float));
          }
    };
}

#endif // __NOMINAL_HH__
